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Abstract. The dynamics of the explosive burning process is highly sensitive to the flame speed model in numerical simulations 
of type la supernovae. Based upon the hypothesis that the effective flame speed is determined by the unresolved turbulent 
velocity fluctuations, we employ a new subgrid scale model which includes a localised treatment of the energy transfer through 
the turbulence cascade in combination with semi-statistical closures for the dissipation and non-local transport of turbulence 
energy. In addition, subgrid scale buoyancy effects are included. In the limit of negligible energy transfer and transport, the 
dynamical model reduces to the Sharp-Wheeler relation. According to our findings, the Sharp-Wheeler relation is insuffcient 
to account for the complicated turbulent dynamics of flames in thermonuclear supernovae. The application of a co-moving grid 
technique enables us to achieve very high spatial resolution in the burning region. Turbulence is produced mostly at the flame 
surface and in the interior ash regions. Consequently, there is a pronounced anisotropy in the vicinity of the flame fronts. The 
localised subgrid scale model predicts significantly enhanced energy generation and less unburnt carbon and oxygen at low 
velocities compared to earlier simulations. 

Key words. Stars: supernovae: general - Hydrodynamics - Turbulence - Convection - Methods: numerical 



1. Introduction 

For supernovae of type la, IHovle & Fowler! {i960) proposed a 
thermonuclear runaway initiated in C+O white dwarfs close 
to the Chandrasekhar limit as the cause of the explosion. 
Since the original proposal, there has been vivid controversy 
of how such an explosion might come about and what the ex- 
act physical mechanism could be. Today the computational 
facilities to process three-dimensional large-eddy simulations 
(LES) of the explosion event are available. Remarkably, these 
powerful means have not aided in arriving at a consensus yet. 
The disagreement stems from some crucial questions. Firstly, 
what is the appropriate flame speed model? Secondly, does 
the explosion completely proceed as a deflagration, or does 
a transition to a delayed detonation set in at some point? 
The d eflagration to detonation transition (DDT) proposed by 
K^okf ikwf ( fl99ll and Wooslev & Weaver! ill 994 appears to re- 
solve the drawbacks of the pure deflagration model. In par- 
ticular, the energy output obtained from simulations with ar- 
tificial DDT is closer to the o bserved one, and le s s car- 
bon and oxygen is left be hind ( Ga mezo et alJ I20041 |2005; 
iGolombek & Niemeverll2005l) . For the theoretical understand- 
ing of thermonuclear supernovae, however, the lack of a con- 
vincing explanation f or the initia t ion of the transition is unsat- 
isfactory dKhokhlov et al.lll997t INiemever & Wooslev! 1 1997t 



!Niemeverll999llZingale et"ail2 005). In the aforementioned nu- 
merical models, a DDT is artifically triggered at more or less 
arbitrary instants of time. 

As for the flame speed model, the controversy is whether 
subgrid scale (SGS) turbulence is mostly driven by Rayleigh- 
Taylor instabilities or dominated by the energy transfer through 
the turbulence cascade. The former point of view holds that the 
magnitude of SGS velocity fluc tuations v' is basically given 
by th e Sharp- Wheeler relation JPavies & Tavlorlll950t [sharp 
1984) 

vrt(0 = 0.5V^ (1) 

where vrt(/) is the asymptotic rise velocity of a perturbation 
of size / due to buoyancy. The effective gravity g e ff is deter- 
mined by the density contrast at the interface between burned 
and unburned material. Setting the turbulent flame speed equal 
to vrt(A), where A is the resolution of the numerical grid, 
has been used in some simulations of type la supernova e 
dGamezo et al.ll2003t ICalder et alJl2003l: iGamezo et al1l2004 . 
How ever, simple scaling argum e nts disfavour this propositio n 
dNiemever & Hillebrandt! 1 19951 INiemever & Kersteinl Il997t) . 
Assuming that non-linear interactions between turbulent ed- 
dies of different size I set up a Kolmogorov spectrum, the 
root mean square turbulent velocity fluctuations obey the scal- 
ing law v'(0 K l 1 ^- Since the Sharp- Wheeler relation implies 



2 



W. Schmidt et al.: A localised subgrid scale model for fluid dynamical simulations II 



Vrt(Z) °c / 1/2 , we have v RT (0/v'(0 I 116 -» towards de- 
creas ing length scales. Consequently. iNiemever & Hillebrandl 
Jl995h adopted a SGS model based on the dynamical equa- 
tion for the turbulence energy fc™, i.e. the kinetic energy of 

IT — — II "3 

unresolved velocity fluctuations (Schumann 1975). The major 
weakness of their approach arises from the fairly tentative clo- 
sures which were formulated ad hoc for LES of stellar convec- 



tion (Clement 1993). 

Various refutations of the scaling argument have been put 
forward. To begin with, the spectrum of turbulence energy 
might be different from the Kolmogorov spectrum. However, 
recent direct numerical simulations support the hypothesis of 
a Kolmo gorov spec trum in buoyancy-driven turbulent combus- 
tion ( Zingale et al 120051) . A more serious concern is that there 
might be not enough time to reach the state of developed tur- 
bulence with a Kolmogorov spectrum in the transient scenario 
of a supernova explosion. This question is difficult to settle a 
priori. For this reason, we took an unbiased point of view and 
accommodated buoyancy effects in the form of an Archimedian 
force term in the SGS turbulence energy model. 

In contrast to the previously used SGS turbulence energy 
model, the new localised model which is thoroughly discussed 
in paper I neither presumes isotropy nor a certain turbulence en- 
ergy spectrum function. This is possible by virtue of a dynam- 
ical procedure for the determination of the SGS eddy-v i scosity 
v sg s = CyAklgs which was adapted from iKim et al] (1999). 
Further more, we apply the co-expanding grid introduced by 
R6pk3 J2005h . The growth of the cutoff length A due to the 
grid expansion poses a challenge for the SGS model because 
of the partitioning between resolved energy and SGS energy 
changes in time. We will show that this rescaling effect can be 
taken into account by utilising the dynamical procedure for the 
calculation of eddy-viscosity parameter C v . The rescaling al- 
gorithm as well as the computation of the Archimedian force is 
explained in Sect. [21 followed by the discussion of results from 
three-dimension numerical simulations in Sect. [3] It is demon- 
strated that the newly proposed SGS model substantially alters 
the predictions of the deflagration model. In particular, we will 
analyse the significance of SGS buoyancy affects. 

2. The flame speed model 

For the relation between the turbulent flame speed s t and the 
SGS turbulence velocity g sgs , we adhere to the results found by 
IPocheaul J 19941) from a theoretical analysis and set 



St = Sla 




(2) 



where C t = 4/3. In the asymptotic regi me of turbulen t burning, 
s t — 2q sgs / V3 which is consistent with lPetersI Jl999l) . 

The evolution of <7 sgs is given by a non-linear partial differ- 
ential equation which is obtained by dividing the equation for 
the specific SGS turbulence energy £ sgs = \q^ g% (see Sect. 3 
of part I) by q sgs . For turbulence driven by the Rayleigh-Taylor 
instability, an additional source term stems from buoyancy ef- 
fects on subgrid scales. This Archimedian force, which is pro- 
portional to the effective gravity due to the density contrast be- 
tween nuclear fuel and ash, directly induces turbulent velocity 



fluctuations. The form of the Archimedian force will be pro- 
posed in Sect. 12. II Moreover, the Eulerian time derivative must 
account for the rescaling of the turbulence energy due to the 
temporal shift of the cutoff length. We will denote the partial 
derivative with respect to the rescaled quantities by d* . The 
Lagrangian time derivative thus becomes 



D* 
Dt 



d*_ 



+ v ■ V. 



(3) 



The complete dynamical equation for the SGS turbulent ve- 
locity is 



D* 1 

-j^-<7s g s--V • (pe K q sgs Vq sgs ) - 4|V? sgs |' 



2 

1 n , n | C *|2 7 j *?SgS 

—c AgeS + e v \s i ~- qsgs d-—. 



(4) 



The rate-of-strain scalar \S*\ is defined by \S*\ 2 = 2S*S* : , 
where S is the trace-free part of the symmetrised Jacobian ma- 
trix of the velocity field, Sy = ^(SyVf + SfVy), and d = S a = <9,v, 
is the divergence. The characteristic length scales € K , l v and € e 
are related to SGS turbulent transport, the rate of energy trans- 
fer from resolved toward subgrid scales and the rate of vis- 
cous dissipation. Each characteristic length can be expressed 
in terms of the effective cutoff length A e ff and a similarity pa- 
rameter: 



C v A e ff 



2 V2A, 



elf 



c f 



(5) 



For the determination of the closure parameters C y , C £ and 
C K see Sect. 4 of part I. The advection of q sgs by the re- 
solved flow is computed with the piece-wise parabolic method 
dColella & Woodwardlll984h . Due to the dissipative effects of 
this numerical sche me on the smallest re solved length scales, 
we set A e ff w 1.6A (Schmidt et al. 2005b). The diffusion terms 
on the left hand side of equation |4] is computed by means of 
fourth order centred differences, and for the source term on 
the right-hand side a semi-implicit Adams-Moulton method is 
used. In the remainder of this section, we describe the calcula- 
tion of the SGS Archimedian force and the rescaling procedure. 

2.1. Archimedian production 

There has been an ongoing debate whether the production of 
SGS turbulence is dominated by the buoyancy of SGS per- 
turbations in the interface between ash and fuel or by eddies 
produced through the turbulence cascade. In the first case, the 
source of energy is the gravitational potential energy, whereas 
non-linear transfer supplies kinetic energy from larger scales 
in the second case. In general, it is quite difficult to sepa- 
rate the energy injection caused by gravity in wave number 
space because gravitational effects are genuinely non-local. For 
Rayleigh-Taylor driven turbulence, however, we know the sim- 
ple Sharp- Wheeler scaling relation Q which provides an al- 
gebraic closure for the Archimedian force density F sgs intro- 
duced in Sect. 3 of part I. Therefore, we propose a novel ap- 
proach which combines both the production of SGS turbulence 
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through the cascade and the Rayleigh-Taylor mechanism in the 
dynamical equation for g sgs . 

Because r sgs has the dimension of an acceleration times 
velocity, we interpret this term as the product of a specific 
Archimedian force and the SGS turbulence velocity q sgs . This 
means the following: In any finite- volume cell those portions 
of the fluid with density less than the smoothed density p will 
experience buoyancy relative to the other portions of higher 
density in the mean gravitational field. In subsonic turbulent 
flows, the small random density fluctuations caused by com- 
pression and rarefaction are expected to produce very little net 
buoyancy. However, a special situation is encountered in the 
cells intersected by the flame fronts. Since the flames are far 
from being completely resolved in numerical simulations of 
SNe la, the substructure of the front in combination with the 
density gradient across the front will produce SGS buoyancy. 
Equivalently, one can think of of perturbations in the flame 
front on length scales Af p < I < A as being Rayleigh-Taylor un- 
stable and producing SGS turbulence. The^re polishing length 
Afp then mark s the lowe r threshold for perturbations to grow 
(see Zingale et al. 2005). Once turbulence has developed, Afp 
can be identified with the Gibson length lc, i.e. the smallest 
length scale on which the flame propagation is affected by tur- 
bulent eddies. Perturbations of size / > A, on the other hand, 
set fluid into motion on numerically resolved length scales. But 
the transfer of energy through the turbulence cascade eventu- 
ally produces SGS turbulence as well. This production channel 
corresponds to the term £ sgs in the equation for the SGS turbu- 
lence energy (see sect 3 of part I). 

The Archimedian force generated by the density gradient 
across flame fronts is given by the effective gravity 



<eff 



Atg, 



where the Atwood number 

At : 



Pf -Pb 



(6) 



(7) 



Pf +Pb 

is a measure for the density contrast between burned material 
and nuclear fuel. In the vicinity of the flame front, the lowest 
order estimate for the SGS buoyancy term is r sgs ~ pq sgs g e ff, 
provided that A ^ Af p . If A becomes smaller than Afp, SGS per- 
turbations in the flame front are not subject to the RT-instability 
and r sgs vanishes. In conclusion, we propose the following clo- 
sure: 

--C A pg e sq sgs - (8) 



r 

1 sgs 



V2 



Here the effective gravity is more precisely defined by 
Seff = X±s(G = O)0(A - A fp )At(p)g, 



(9) 



where x±nh{G = 0) is the characteristic function of all cells 
for which the distance from the flame front (represented by 
G(x, t) = 0) is less than 6, and 8 is the Heaviside step func- 



tion, i.e. 0(A - Afp) 



1 for A > Afp and zero otherwise. 



The 



Atwood number is expressed as a function of the mean den- 
sity which is obtained fro m a fit to the numerical data from 
Tim mes & Wooslevl ( 1 19921) : 



At(p) = \ 



0.0522 + 



0.145 0.0100 



P9 



(10) 



The closure © does not include all of the intricate effects 
that contribute to SGS buoyancy. It merely captures what is 
presumably the leading order effect. In fact, one would have to 
model the interaction between turbulent potential and kinetic 
energy fluctuations on unresolved scales. Unfortunately, there 
exists no theoretical framework for this task yet. Moreover, the 
concept of a SGS Archimedian force entails a violation of en- 
ergy conservation because the contribution of r sgs to the pro- 
duction of turbulence energy is not balanced in the resolved 
energy budget. Consequently, the total energy of the system ef- 
fectively increases. However, r sgs is non-zero only in a small 
volume fraction. For this reason, the resulting violation of en- 
ergy remains negligible relative to the total energy budget. 

In order to determine the parameter Ca, we observe that 
the Sharp- Wheeler SGS model is obtained as an asymptotic 
relation in the limiting case of neglecting the non-local trans- 
port I) sgs , the turbulent energy transfer £ sgs and the pressure- 
dilatation /l sgs (see Sect. 2 of paper I for definitions). Dropping 
the corresponding terms in equation one obtains 



d i _ 



(11) 



for a fluid parcel in the vicinity of the flame front. In the sta- 
tionary regime, this equation has the fixed point solution: 



2CAA e ffg e ff , . . 

<7sgs - -y ^ - VRT(A e ff,). 



(12) 



Consistency with the Sharp-Wheeler relation Q implies Ca = 
C f /8. Since C e » .5 ... 1 .0 for developed turbulence (see 
ISchmidt et all2005al) . the estimate Ca ~ 0.1 is obtained. 

2.2. Rescaling of the subgrid scale turbulence energy 

If a non-static, co-moving grid is used, the implicit filter 
( ) e ff introduced in paper I, Sect. 3, becomes time-dependent. 
Therefore, time-derivates do not commute with the operation of 
filtering and additional terms arise in the dynamical equations. 
These terms are equivalent to the additional fluxes which are 
included i n the imp lementati on of the Riemann solver for mov- 
ing grids llMullerl 1 994l lRonkel2005l) . However, there is a sub- 
tlety related to the shifting cutoff which separates the resolved 
and unresolved scales. As the grid expands homologously with 
the bulk of the white dwarf, the grid resolution A gradually 
decreases in time and the growing cutoff length entails a grad- 
ual rise of the SGS turbulence energy. This rise is inherent to 
the grid geometry and immediately affects the decomposition 
of the energy budget. The two- thirds law fo r developed turbu- 
lence implies (<7 sgs ) °c A I//3 (see Frisch 1995] Sect. 5). Thus, it 
is easy to rescale the mean value of q sgs if A changes by a small 
fraction 6 A/ A: 



1 6A\ 



(13) 



Because this scaling law is a statistical rule one cannot ex- 
pect that it holds locally. However, the dynamical procedure for 
the calculation of the eddy-viscosity also can be utilised for a 
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localised rescaling law. Let A, be the grid resolution at time t, 
and Af_is ( the slightly smaller resolution prior to the last time 
step. Applying the test filter introduced in paper I, Sect. 4.1, at 
time t, the turbulent velocity q T associated with velocity fluctu- 
ations on length scales greater than /3A, and smaller than yjflA, 
is obtained. Here /3 - A e fr /A is the constant ratio of the effec- 
tive cutoff length to the size of the grid cells. Since the Riemann 
solver does not account for the fractional growth of the turbu- 
lence energy due to the shift of the cutoff, the result of advanc- 
ing the dynamical equations from t — St to t is q^' ' . Now an 
estimate of qf^ s can be made locally via interpolation of the 
turbulence energy in length scale space. Using the contracted 
Germano identity (paper I, Sect. 4.1), the total turbulence en- 
ergy associated with the test filter length at time t is given by 



turb 



(14) 



The unknown is the rescaled SGS turbulence energy k S g S . 
Neglecting compressibility effects and variations on sub-test 
filter lengths, we set k^' - kfg S + k m . Then linear interpo- 
lation between A,_,5, and yxA, yields: 



iA ~ tJk-st i t u 

ft sgs — "sgs ^ j _ y T ' 

where the interpolating factor / is given by 



A, - A,s, 
y T A, - A,_ 



(15) 



(16) 



Due to the smallness of a CFL time step, the fractional changes 
of the cutoff length will be small. Hence, / «: 1. 

The problem with the rescaling law J15l > is that it fails to ac- 
count for the correct asymptotic behaviour in the limit of fully 
developed turbulence. On account of the Germano identity (pa- 
per I, Sect. 4.1), one would expect the statistical relation 



!)<*«.> = <£ (T) > 



(17) 



for regions of nearly homogeneous turbulence obeying 
Kolmogorov scaling. This relation is asymptotically repro- 
duced by the rescaling modified law 



rfd-/) 



A-si 



f 



2/3 



■/ 



2/3 



f 



(18) 



which results from the interpolation of the turbulence energy 
divided by the associated length to the power 2/3. This is the 
rescaling law which was implemented for the numerical simu- 
lations of thermonuclear supernova explosions discussed in the 
next section. 

3. Numerical simulations 

In the following, we will present results from several 
numerical simulation s usin g the methodology outlined in 
Ropke & Hille brandtl (2005J). In essence, the piece-wise 
parab olic method is used to solve the hyd rodynamical equa- 
tions (Colella & Woodward 1984: iFrvxell et alJfl989l) and the 
evolution of the flame-fronts is computed by means of the level 



set method in the passive i mplementation dOsher & Se thian 
1988| iReinecke et alJ ^999). The implemented e quation of 
state f or electron-degenerate matter is described in IReinecke] 
(12001 1) . Sect. 3.2. Thermonuc lear burning is m odelled by sim- 
ple representative reactions dReinecke et al.l 2002): 12 C and 
16 C is fused to 56 Ni and a-particles at densities higher than 
5.25 • 10 7 g cirT 3 , whereas 24 Mg is produced at lower densities 
in the late stage of the explosion. Finally, all reactions cease 
below 10 7 gcirT 3 . This threshold presumably marks the tran- 
sition from the flamel et to the broken reaction zo ne regime of 
turbulent burning (see Nieme ver & Kersteinl 1997b . The correct 
treatment of the burning process in this regime is still a matter 
of debate and, for this reason, not included in the present sim- 
ulations. 

As initial model, we choose a white dwarf of mass M = 
1 AMq composed of equal mass fractions of carbon and oxy- 
gen with a central density p c = 2.0 • 10 9 g cm' 3 and temperatur e 
T c = 7.55 ■ 10 8 K. As suggested bv lWunsch & WooslevN2004 . 
the radial temperature profile is given by a parabola with a cut- 
off at the thermal radius A = 7.35 • 10 7 cm: 



T(r) = T c 



-(xi 



9(r - A) + T o 0(A - r), (19) 



where 9 denotes the Heaviside step function. The thermal ra- 
dius specifies the size of the convective core prior to the run- 
away. At larger radii, the matter is isothermal with To — 5 • 
10 5 K. In the centre, we set an axisymm etric initial bur ning 
region with sinusoidal perturbations (see lRopke & Hillebrandtl 
2005). In order to achieve higher resolution, only single oc- 
tants subject to reflecting boundary conditions were evolved in 
the simulations discussed here. 

Moreover , we applied the co-expanding grid technique of 
lR6pkd J2005). Thereby, it is possible to maintain an equidistant 
grid geometry over the whole domain of turbulent burning at 
any stage of the explosion, even when the ejecta have expanded 
by a large factor compa red to the initial size of the white dwarf. 
Recently, iRopke et al.l d2005t) have combined this t echnique 
with the grid geometry used by R einecke et alJ J2002I) in order 
to capture the burning process in the interior with optimal res- 
olution, while using a coarser grid with exponentially increas- 
ing cells outside. The hybrid grid, even at moderate resolution, 
enables us to resolve details of the turbulent flame dynamics 
which used to be inaccessible for non-adaptive schemes. All 
numerical simulations presented in this article feature a hybrid 
grid. 

The non-uniform grid geometry poses certain difficulties 
when applying the localised SGS model. In Sect. 12.21 we 
showed that it is relatively easy to account for the variation in 
the time domain due the co-expansion of the grid. In the case of 
non-uniform grids, however, the filter operation does not com- 
mute with spatial derivatives in the dynamical equations. Apart 
from that, the weighing of nodes for the discrete filtering pro- 
cedure becomes dependent on the location. This would lead 
to substantial complications in the numerical implementation. 
Fortunately, we found a simple solution: Since the turbulent 
dynamics mostly takes place in the burning region which is 
contained within the uniform part of the grid, we computed the 
eddy-viscosity parameter C v only in this region and set the rate 
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Fig. 1. Time evolution of integrated quantities for three simulations with identical initial conditions and resolution 256 . In one 
case, Clement's SGS model with wall proximity functions (WPF) was used. For the other two simulations we applied the localised 
SGS model with Archimedian production, respectively, switched off and on. 



of turbulent energy transfer equal to zero in the exterior. Later 
in this section, we will demonstrate that neglecting the energy 
transfer outside the burning region can be justified a posteriori. 

The original SGS turb ulence energy model implemented by 
iNiemever & H illebrand 3 \ 19951) is based on statistical closure 
parameters. lClemenll i ll 9931) suggested to set C v — 0. 1 W, where 
W is an empirical wall proximity function (WPF), 



W = min 



100,max|0.1,-l(T 4 ^ 



(20) 



Since e ml ~ cj, the ratio e- m lqi K% is basically an inverse Mach 
number squared. If there is little turbulence energy, C v is con- 
siderably enhanced. On the other hand, C v becomes smaller 
than 0.1 if the SGS turbulence velocity exceeds a few percent 
of the speed of sound. This behaviour of the eddy-viscosity 
parameter is qualitatively different from the prediction of the 



dynamical procedure which implies less energy transfer if tur- 
bulence is still developing but enhan ced transfer in the fully 
developed case dSchmidt et al.l2005al) . 

This is clearly reflected in the time evolution of the turbu- 
lence energy in two simulations which differ only in the SGS 
model. The graphs of the mass-integrated SGS turbulence en- 
ergy are plotted in the top panel on the left of Fig. [2 There 
are two variants of the localised SGS model, one including 
Archimedian production as described in Sect. 12.11 whereas it 
is assumed that energy transfer through the cascade is the only 
source of turbulence production in the alternative model. In 
contrast to Clement's model, there is initially very little SGS 
turbulence energy followed by a much steeper rise in the simu- 
lations with the localised SGS model. The rapid growth of tur- 
bulence energy can be attributed to the substantially stronger 
turbulence production within the time interval from 0.3 to 1 .0 s 
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Fig. 2. Time evolution of integrated SGS quantities for a series of simulations with varying resolution. 



(see right bottom panel in Fig.[Q. In the second half of the com- 
bustion phase, the total SGS turbulence energy is almost one 
order of a magnitude larger which enhances the flame propa- 
gation speed accordingly. At later times, the discrepancy be- 
comes even more pronounced because the rescaling of £ sgs im- 
plemented in the localised model feeds kinetic energy into the 
subgrid scales against the action of SGS dissipation. The net re- 
sult is a significant enhancement of the explosion energy and a 
larger yield of burning products (see bottom panel on the left of 
Fig-Q- Also note that the additional production of SGS turbu- 
lence by the Archimedian force term in equation @ increases 
the explosion energy even further. Compared to the reference 
simulation with Clement's model, the final kinetic energy of 
0.472 ■ 10 51 erg is about 25 % greater. However, this does not 
imply that Archimedian production dominates over the turbu- 
lence cascade. In fact, the evolution of the SGS turbulence en- 
ergy differs only little, as one can see from the plot in Fig.^ 
In conclusion, SGS buoyancy effects appear to influence the 



explosion but turbulent energy transfer from resolved scales is 
nevertheless the primary source of SGS turbulence production. 

According to the scaling argument mentioned in the intro- 
duction, buoyancy should become even less important relative 
to the turbulent energy transfer with increasing resolution. This 
is indeed observed in a series of simulations with the reso- 
lution varying between 128 3 up to 384 3 grid cells. The time 
evolution of the integrated rate of energy transfer and specific 
Archimedian force, respectively, is plotted in the top panels 
of Fig. |2] In order to interpret these graphs, it is important to 
note that the SGS energy transfer is expected to become sta- 
tistically scale invariant in the case of Kolmogorov turbulence. 
This follows from the scaling law v'(l) oc I 1 / 3 for the turbu- 
lent velocity fluctuations. Hence, k sss oc A 2 ^ 3 . Since the char- 
acteristic time scale of turbulent velocity fluctuations scales 
with Z 2 ' 3 , it follows that the average time derivative of k sgs is 
scale invariant. The rate of energy transfer per unit mass, on 

1/2 9 

the other hand, is proportional to A£ S g S \S*\ . This expression is 
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also scale-invariant because the rate-of-strain scalar \S*\ mea- 
sures the inverse time scale of the smallest resolved velocity 
fluctuations, i.e. |5"*| 2 oc A~ 4 ^ 3 , while Afcjg 2 oc A 4/3 . The com- 
puted rate of energy transfer plotted in the top panel on the left 
of Fig.|2exhibit peak values which are within the same order of 
magnitude, although the energy transfer seems to be underes- 
timated for the lowest resolutions. Initially, the turbulence en- 
ergy rises exponentially at a rate which changes only little with 
resolution (right bottom panel in Fig.|2}. This is reflected in the 
nearly coinciding graphs of the rate of SGS energy transfer up 
to t w 0.3 s shown in the left panel on the top of Fig.|5] 

The SGS turbulence energy in the regime of fully devel- 
oped turbulence decreases for higher resolution. This trend can 
be discerned particularly in the post-burning phase, in which 
no further energy is injected and the turbulent flow begins to 
decay. From the initial production phase to the post-burning 
phase, however, a rather complicated behaviour becomes man- 
ifest as the result of the interplay between turbulence produc- 
tion by the strain of the resolved flow, the SGS Archimedian 
force and the grid expansion. The explosion energetics plotted 
in Fig.|5]shows the following behaviour depending on the nu- 
merical resolution (also see Table [Q: Whereas the final value 
of the total energy is about the same for the lower resolutions, 
there is a significantly enhanced yield of energy in the case 
N = 384 3 . However, this does not imply that the model fails 
to converge with increasing resolution. Turbulent flow regions 
are confined in a fraction of the numerical grid, whereas the 
greater part of the grid is overhead required for modelling the 
non-turbulent outer parts of the expanding star and some por- 
tion of the surrounding quasi-vacuum. In the case N = 256 3 , 
for example, significant SGS turbulence production occurs in 
the inner 100 3 cells at t — 0.5 s. This is the time of maximal in- 
tegrated turbulence production. However, in paper I we demon- 
strated that 100 cells in each spatial dimension is definitely not 
sufficient to resolve developed turbulent flow sufficiently far 
down toward the inertial subrange using PPM. Although this 
is merely a crude estimate, it appears plausible that even the 
supernova simulation with N = 384 3 resolves the turbulent dy- 
namics only marginally. As a further indication, the plateau- 
like flattening of the rate of production and dissipation, respec- 
tively, can be seen in the left panels of Fig. 13 for the highest 
resolution only. We interpret the flattening as a consequence 
of local statistical equilibrium between production and dissipa- 
tion. 

Unfortunately, we were not able to perform a run of still 
higher resolution due to the limitations of our computational 
resources. In any case, we expect that N = 512 3 grid cells in 
one octant would be sufficient for an accurate modelling of the 
turbulent dynamics, whereas simulations with less resolution 
can be utilised to discern t rends in para meter studies. A similar 
conclusion was drawn bv iRopkel 120051) from a series of two- 
dimensional simulations, in which a pronounced jump of the 
total energy was found between the N = 256 2 and the 5 12 2 run, 
respectively, while more or less the same energy was obtained 
for N > 512 2 . Moreover, snapshots of the zero level set for 
varying resolution suggested that secondary Kelvin-Helmholtz 
instabilities are barely or not at all resolved with N < 256 2 . 




0.1 1.0 
t [s] 

Fig. 3. Time evolution of the total energy for the same simula- 
tions as in Fig. 13 

Table 1. Total release of nuclear and kinetic energy and total 
masses of iron group ("Ni") and intermediate mass ("Mg") el- 
ements corresponding to Fig. [3] 



N 


E auc [10 51 erg] 


E km [10 51 erg] 


M Ni /M 


M Mg /M Q 


128 3 


0.963 


0.433 


0.523 


0.175 


192 3 


0.970 


0.442 


0.529 


0.172 


256 3 


1.000 


0.472 


0.548 


0.172 


384 3 


1.087 


0.560 


0.586 


0.206 



Details of the SGS dynamics are illustrated by contour plots 
of two-dimensional spatial sections from the simulation with 
N = 384 3 grid cells (Figs.|H|5] and|6}. In each Fig., the follow- 
ing dynamical terms of equation (@) are plotted: 

1 . Rate of production caused by strain, ( V \S *| 2 (left top panel). 

2. Specific Archimedian force 0. lg e ff (right top panel). 

3. Rate of dissipation -^q sgs d - q\ g Ji e (left bottom panel). 

4. Rate of diffusion ± V ■ (p( K q sgs Vq sgs ) - 4|Vg sgs | 2 (right bot- 
tom panel). 

The flame surface as given by the zero level set is indicated 
by the contours in white. Note that these quantities have the 
dimension of acceleration. Fig. @] shows the typical Rayleigh- 
Taylor mushroom shapes which have formed out of the initial 
sinusoidal perturbations at time t = 0.3 s. Significant energy 
transfer is concentrated in small regions and there is little dis- 
sipation yet. Comparing to Fig. 12 one can see that turbulence 
production is just about to rise. At t — 0.45 s, the rate of energy 
transfer has reached its maximum and is spread all over the in- 
terior of the flames (see Fig.|3}. The acceleration of SGS fluid 
parcels subject to the largest strain exceeds 10 6 times the gravi- 
tational acceleration on Earth relative to the resolved flow. The 
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Fig. 4. Contour sections showing the contributions to the evolution of the SGS turbulence velocity g sgs given by equation 0) at 
t = 0.3 s. Only the inner region of the grid with N = 384 3 cells is shown. The white contours represent the sections through the 
flame surface. 



SGS buoyancy is typically by an order of a magnitude smaller. 
Both dissipation and transport due to SGS turbulent diffusion 
are comparable to the rate of energy transfer at this time. In 
the unburned material outside, on the other hand, there is virtu- 
ally no SGS turbulence. Thereby, it is confirmed that switching 
off the energy transfer terms in the non-uniform gird regions at 
sufficient distance from the flame fronts is a reasonable simpli- 
fication. Obviously, the flow is highly anisotropic in the vicin- 
ity of the flames which highlights the necessity of a localised 
SGS model. In Fig. [6] one can see that turbulent energy transfer 
is declining and becoming small relative to the Archimedian 
force near the flame front. However, this does not imply that 
the amount of SGS turbulence and, thus, the turbulent flame 
speed is dominated by the SGS buoyancy because the bulk of 
SGS turbulence energy has been produced by transfer of ki- 



netic energy through the turbulence cascade and diffusion acts 
to redistribute this energy into regions with little production. 

Fig. 1101 shows the evolution of the SGS turbulence veloc- 
ity q sgs at the flame fronts in three-dimensional visualisations. 
The grid lines roughly indicate the uniform part of the numer- 
ical grid. The corresponding absolute scale is indicated by the 
size X un [ and the corresponding number of cells iV U ni- In the first 
three snapshots one can see the growth of the initial perturba- 
tions. The axial symmetry is gradually broken by the forma- 
tion of secondary instabilities. At t — 0.45 s the smaller plumes 
originating from these instabilities are highly turbulent. From 
t = 0.6 s onwards, the system increasingly looses its memory 
of the initial condition and the turbulence intensity at the flame 
surface is abating and levelling. The last snapshot at t — 1.5 s 
shows a complex structure with features over a wide range of 
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Fig. 5. The same plot as in Fig.|4]at t = 0.45 s. 



scales. There appear to be five or six major modes which even- 
tually prevail. The resulting layering of nuclear species is illus- 
trated by the contour plots of the corresponding mass densities 
in Fig. [7] Nickel is concentrated in the central region, whereas 
both magnesium and unburned carbon and oxygen are found 
further outside. The outermost layers and the narrow down- 
drafts between the convective fingers of nuclear ash are com- 
posed almost exclusively of carbon and oxygen. The stratifica- 
tion of the nuclear species in the explosion ejecta is reflected 
in the corresponding mass density functions dM/dvv, where v r 
is the radial velocity component. In particular, Fig. [8] shows 
that little carbon and oxygen is found for velocities less than 
3000 km s _1 in the late phase of almost homologous expansion. 

To understand the flame dynamics, it is instructive to con- 
sider the probability density function (PDF) of the logarithm of 
the flame propagation speed s t over the surface of the flame. 
The PDFs for several instants of time are plotted in Fig. [9] 



Note that integrating each PDF over the decade logarithm of 
the speed yields unity. Also shown are the PDFs of q sgs and 
VRx(A e ff )■ The relation between the turbulent flame speed s t and 
the SGS turbulence velocity q sgs is formulated in equation (0. 
During the first tenth of a second, s t is basically given by the 
laminar flame speed. Then the flame propagation becomes in- 
creasingly affected by SGS turbulence. From about 0.3 s on- 
wards, s t is dominated by q sgs . At later times, one can see the 
asymptotic relation s t 2q sgs / y3. The Rayleigh-Taylor ve- 
locity scale VRx(A e ff) is initially much smaller than the laminar 
burning velocity. As the flame propagates outwards, both the 
gravity and the density contrast at the flame surface become 
larger and VRj(A e ff) increases. Eventually, the PDF of VRr(A e ff) 
tends toward a rather narrow peak around 10 7 cm s _1 . The PDF 
of q sgs , on the other hand, extends over a substantially wider 
range. For this reason, the localised SGS model generates more 
variation in the propagation of the flame front in comparison to 
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Fig. 6. The same plot as in Fig.|4]at t = 0.75 s. 



the Sharp-Wheeler model. This is expected because the Sharp- 
Wheeler relation is ignorant of the interaction between sub- 
grid and resolved scales and the effects of non-local transport. 
Nevertheless, VRT(A e ff) is seemingly a velocity scale which is 
representative for the magnitude of q sgs at the flame surface 
during most of the burning process. 

4. Conclusion 

We applied the SGS turbulence energy model to the large eddy 
simulation of turbulent deflagration in thermonuclear super- 
nova explosions. The novel features of this model are a lo- 
calised closure for the rate of energy transfer, an additional 
Archimedian force term which accounts for buoyancy effects 
on unresolved scales and the rescaling of the SGS turbulence 
energy due to the shift of the cutoff length in simulations with 
a co-expanding grid. We found that the production of turbu- 
lence is largely confined to the regions near the flame fronts and 



in the interior ash regions. Consequently, there is pronounced 
anisotropy at the flame surface which can be tackled by the 
localised SGS model only. The Archimedian force contributes 
noticeably to the turbulent flame speed, particularly once the 
flame surface has grown substantially. However, the dominat- 
ing effect is the energy transfer through the turbulence cascade. 
In the late stage of the explosion, sustained turbulence energy 
comes from the rescaling, while the major dynamical contribu- 
tion is SGS dissipation. Furthermore, it appears that numerical 
grids with more than N = 256 3 cells in one octant are neces- 
sary in order to sufficiently resolve the turbulent dynamics in 
the burning regions and to obtain converged results. 

An investigation of probability density functions over the 
flame fronts (see Fig. |9} reveals that the Rayleigh-Taylor ve- 
locity scale v'R-r(A e jf) given by the Sharp- Wheeler relation Q 
is not negligible compared to the SGS turbulence velocity 
q sgs , once the regime of fully turbulent burning has been en- 
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Fig. 7. Total and partial mass densities for the same simulation as in at Fig. gjjfjat t = 1.5 s. 



tered. This reflects the slow decrease of the ratio VRT(A e (f)/q , sgs 
with the numerical cutoff scale according to the scaling argu- 
ment discussed in the introduction. The underlying scaling re- 
lations do not necessarily apply to transient and inhomogenous 
flows as in the supernova explosion scenario. The PDF for q sgs 
shows a considerably wider spread than the sharply peaked 
PDF for VR-r(A e jf). We interpret this observation as a conse- 
quence of the additional physics in the localised SGS model, 
which also encompasses turbulent energy transfer (i.e. interac- 
tion between resolved and subgrid scales) and turbulent trans- 
port (i.e. non-local interactions among subgrid scales). The re- 
lation between the Sharp- Wheeler and the turbulence energy 
models may be analogous to the relation between the mixing 
length and Reynolds stress models of convection. 

The final kinetic energy in the simulation with the high- 
est resolution is about 6 • 10 50 erg. The produced mass of 
iron group elements, O.58M0, falls within the range deduced 



from observations of type la supernovae jLeibundgud EoOO). 
However, some observed events are substantially more ener- 
getic. Regarding the numerical simulations, the explosion en- 
ergy is very sensitive to the initial conditions and the localised 
SGS model appears to increase the sensitivity even further. 
Using initial conditions that are different to the highly artificial 
centrally ignited flame is in progress. A persistent difficulty is 
the large amount of left over carbon and oxygen (see Fig. Q 
and[S}. It is not clear yet to what extent this problem can be 
solved with the aid of the localised SGS model in simulations 
with more realistic ignition scenarios. It appears more likely 
that a different mode of burning is required in the late explo- 
sion phase. The currently implemented numerical burning ex- 
tinction at a density threshold of 10 7 g cm -3 is mostly arbitrary 
and should be replaced by a physical criterion motivated by the 
properties of distributed burning at low densities. This might 
turn out to be an alternative to the DDT scenario. 
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Fig. 8. Density functions of mass in radial velocity space for 
the major nuclear species in the same simulation as in Fig.|4j{7] 
at t = 5.0 s. 
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Fig. 9. Probability density functions of VR-r(A e ff), <? S gs and s t 
over the flame surface at several instants of time. 
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fig.c3_3d_015.png 
(a) t = 0.15 s, N um = 225 3 , X um = 320km 



fig_c3_3d_030.png 

(b) t = 0.3 s, N um = 242 3 , X um = 484km 



fig_c3_3d_045 .png 

(c) t = 0.45 s, iV uni = 258 3 , X um = 721 km 



fig_c3_3d_060.png 

(d) t = 0.6 s, 2Vu„i = 283 3 , X um = 1260 km 



fig_c3_3d_075 .png fig_c3_3d_150.png 

(e) t = 0.75 s, yV um = 308 3 , X um = 2190 km (f) r = 1.5 s, = 340 3 , X um = 9120 km 

Fig. 10. Evolution of the flames in the simulation with N = 384 3 grid cells. The colour shading indicates the value of g sgs on a 
logarithmic scale. 



